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Section 1 


INTRODUCTION 

A structural engineer is often faced with the need to predict the 
response of structures under failure conditions in attempting to make 
them fail-safe or crashworthy. Under failure conditions the load magni- 
tudes may be large and the structural response may involve large perma- 
nent deformations, rupture and tearing. 

With finite deformations the equations of motion of the structure 
are coupled in products of the derivatives of displacements and stresses: 
the unknowns of the problem. Material yielding may introduce non-linear 
stress-strain relations. Thus, nonlinear transient response prediction 
requires the evaluation of deflections, velocities, accelerations, 
stresses and strains in a structure of complex geometry and ductile 
material under time-varying loads that may cause the members of the 
structure to undergo large motions and deformations and/or respond 
plastically . 

A number of well-known computer codes which can be used for post- 
impact studies of aircraft may be identified [1,2], They can be classi- 
fied as general or special purpose programs. General purpose programs 
such as MARC and NASTRAN^ include some of the capabilities necessary for 
crash simulations. Special purpose programs such as the present computer 
program ACTION, provide in addition special modeling and evaluation 
processes which are particularly useful for crash. Because special purpose 
programs focus on a particular problem class, they are potentially more 
efficient for that class, and presumably require less computer resources, 
and are more manageable. 

In the partitioned spectrum of crash simulators defined by Mclvor, 

^NASTRAN: Registered trademark of the National Aeronautics and Space Administration. 



[1], ACTION qualifies as a level four code: providing the capability to 

represent the material and geometric nonlinear transient behavior of a 
structure composed of truss, frame and membrane elements. In com- 
parison with other level four simulators, ACTION offers the following 
unique features: 

1. The ability to automatically control time discretization error. 

2. Representation of impact with a rigid barrier including treat- 
ment of gapping and friction effects. 

3. Logic specifically designed to minimize data transfers by keep- 
ing all working data in core and providing an analysis mode which avoids 
generation of large stiffness matrices in explicit form. 

4. Response data defining the allocation of stored and dissipated 
energies . 

The ACTION code, although quite suitable for nonlinear static 
analysis is primarily designed for analyzing response of vehicles crash- 
ing into a rigid or a deformable barrier. The objective is to predict 
analytically the response of a lightweight aircraft during a crash. The 
simulation is based on a discretized model of the structure using finite 
elements which may undergo finite motions and deformations and may 
respond plastically. A transient analysis of such a model then yields 
the displacements, velocities, accelerations, internal loads and 
stresses, at points of interest, under time varying loads that may cause 
complete failure of the structure. 

The objective of such a nonlinear transient analysis is to develop 
an understanding of the multi-faceted relationship between the complex 
structural configuration of an aircraft and its response during crash. 
Such an understanding can provide the basis for crashworthy design of 
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lightweight aircraft, better restraint systems and efficient energy 
absorption devices which will reduce passenger trauma. 

This report outlines the theoretical basis for the ACTION computer 
code. Instructions for the preparation of input and the interpretation 
of results can be found in reference [3] . 
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Section 2 


FORMULATION OF ACTION 

This section describes the basis and highlights the intricacies of 
nonlinear transient analysis of a structure. Because of the complexity 
of structural geometry and response it is necessary to construct an 
approximate mathematical model. The objective is to predict the dis- 
placements, velocities, accelerations, internal loads and stresses at 
points of interest under time-varying loads that may cause complete 
failure of the structural model. 

The mathematical model is a finite element displacement model. 
Simulation consists of discretizing the actual structure by finite ele- 
ments, approximating the response of each element by a finite number of 
deformation states expressed as linear functions of generalized joint 
displacements and analyzing the mathematical model numerically. 

Most researchers prefer to use a step-by-step incrementally linear 
approach for the solution of the nonlinear equations of the finite ele- 
ment model. In this process, sometimes referred to as the vector ap- 
proach, the load is applied in increments which are sufficiently small 
so that a linear analysis will approximately represent the structural 
response for in increment. An iterative technique at constant load is 
used to satisfy equilibrium exactly for the partial load. The deformed 
geometry and stress state for a partial load is used as the initial 
state for the next load increment. Turner et al. [4], in their pioneer- 
ing paper, provide a description of the process for finite deflection 
analysis by the step-by-step technique. Martin [5] provides some illus- 
trations. Oden [6] and Hibbitt et al. [7] show applications to problems 
involving large strains and finite displacements. A comprehensive re- 
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view of the many abberations of the incremental method for geometrically 
nonlinear problems is provided by Haisler et al [8]. Armen et al. [9], 
[10] Belytschko [11], [12] generalize the step-by-step process to pro- 
blems involving both geometric and material nonlinearities. The work 
of these authors represents the many applications of the step-by-step 
approach that may be found in the literature. 

ACTION represents the structural characteristics by a scalar func- 
tion. The equilibrium configuration is established based on derivatives 
of the function. This approach, sometimes referred to as the scalar ap- 
proach, has been successfully used in nonlinear structural analysis by 
Bogner et al, [13]. Mallet and Berke [14] show results of applying the 
process to truss structures. For transient analysis the formulation 
couples a step-by-step numerical integration in time with the function 
minimization procedure. Discrete steps are taken in the time co-ordinate 
and function minimization is used at each discrete time point to find 
the solution. Young [15] was probably the first to illustrate the use 
of this approach for nonlinear transient response simulation. 

The approach consists briefly of the following steps: 

1. Assume suitable element displacement field as a function of 
the local spatial co-ordinates with generalized joint dis- 
placements as unknown coefficients. 

2. Relate the element generalized displacements to global gen- 
eralized displacements of the system accounting for pre- 
scribed boundary conditions. 

3. Develop expressions for strains in each element as functions 
of the global generalized displacements. (These expressions 
will account for nonlinearities in the derivatives of the dis- 
placement field.) 
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Determine the corresponding stresses and strain energy den- 
sities in each of the elements of the assemblage using respec- 
tive material models. 

5. Integrate the strain energy density over the volume of each 
element to yield element strain energy as a function of the 
global generalized displacements. 

6. For transient analysis, knowing the values of the generalized 
displacements, velocities, accelerations, etc., at time t 
extrapolate assuming displacements vary as prescribed func- 
tions of time. 

7. Determine by energy minimization, the configuration which im- 
plies satisfaction of equilibrium at the end of the next (load) 
time increment. 

The simulation provides accurate results for static analysis with 
linear material behavior. For transient analysis with inelastic material 
behavior several factors limit the accuracy: 

1. Because of the assumed displacement-time relations for the 
response the approximation does not accurately describe the 
response for large time steps. Depending upon the type of the 
temporal integration scheme used, there is a possibility of 
divergence of the solution from the true solution as time steps 
accumulate. 

2. For inelastic stresses, the evaluation of the strain energy 
may be in error for certain deformation states due to approxi- 
mations in integration. This error may be reduced by modeling 
the affected element with a number of smaller elements. 

3. Because of discretization of the actual structure by an assemb- 
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lage of kinematically admissible finite elements, in general, 
the calculated solution can be expected to be in error from 
the true solution for the actual structure due to large geo- 
metry changes of particular elements. This error can also be 
reduced by modeling the given structure with more elements. 

4. The inelastic material representation may be inadequate in some 
cases because the constitutive equations of ACTION are rather 
simple idealizations of true behavior. 


7 



Section 3 


DEFORMATION MODEL 

The deformation model is a mathematical model which characterizes 
the deformation of the structure in terms of unknown variables. The 
model is synthesized from deformation states of each element of the 
structure. These deformations are expressed in terms of displacements 
of "joints" of the system, points at which the elements interface. 

The displacement field within each element is chosen as a continu- 
ously differentiable function of the local spatial co-ordinates and the 
joint displacements. The field maintains interelement continuity of the 
essential derivatives and includes constant strain states so that the 
representation provides a Galerkin model of the system. The generalized 
joint displacements of each element are related to the global displace- 
ments of the assemblage of the elements. These relations, which can be 
interpreted as transformations of the local generalized displacements, 
may be linear or nonlinear depending upon whether the motions and defor- 
mations of the elements are infinitesimally small or finite. For large 
angular changes, these transformations are accomplished using Euler 
angles which are linearly independent by virtue of the fact that the 
rotations are performed in a prescribed order. 

The remainder of this section formulates the deformation character- 
istics of each type of element used in the ACTION simulator and the 
transformations to relate element behavior to the common rectangular 
cartesian coordinate system. 

3.1 TRUSS ELEMENT [15] 

A truss or a rod element is a structural component of uniform cross- 
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section which is initially straight. It is assumed that the element is 
capable of resisting only axial deformation imposed by joint displace- 
ments through momentless connections. Using a linear interpolation 
function for the axial displacement, closed form expressions for stresses 
and strain energy for the linear elastic range are developed. For the 
inelastic range development of the expression for the strain energy 
involves the use of the material model of the element. 

3.1.1 Deformation 

Figure 3-1 shows the initial and deformed positions of a typical 
truss element specified by the co-ordinates of the two joints, p and q, 
which it connects. These co-ordinates are defined with respect to a 
fixed global system of reference axes X, Y and Z. x and y denote the 
local co-ordinates axes corresponding to the deformed nodal configura- 
tion. Motion of the element is defined by finite displacement vectors 
Up and U^ of the two joints. The components of the displacement vector 
IJ in the X, Y, Z co-ordinate directions are denoted by U, V and W re- 
spectively. Element deformations are denoted by an axial displacement 
function u(x) . It is convenient to express the element strains with 
reference to the deformed nodal configuration, finite nodal displace- 
ments being accounted for in the transformation from global to the local 
displacements , 

From Fig. 3-1 the initial length, L, of the element is given by 
L = [ (X - X ) 2 4- (Y - Y ) 2 + (Z - Z ) 2 ] 1/2 (3-1) 

q p ^ q p q P 

The deformed length, L, is given by 

L = [ (X + U -X -U) 2 +(Y + V -Y -V) 2 +(Z +W -Z -W ) 2 ] 1/2 (3-2) 

qqpp qqpp qqpp 

Hence, the change in length, DL, is given by 
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DL = L[1 + 


2 (AXAU + AYAV + AZAW) 


AU 2 + 


AV + AW ,1/2 
o J 


- L (3-3) 


where A is the difference operator for q and p end values. 

Equation (3-3) gives the change in length for joint displacements 
of any magnitude. The quantity inside the brackets in Eq. (3-3) should 
be positive but due to manipulation errors, values less than zero may be 
realized. Then DL is assumed to be -L: the maximum value physically 

possible. For values of the quantity inside the brackets close to one, 
DL is evaluated by using a binomial expansion. 

Next, the assumption of the usual linear interpolation function in 
the corotational co-ordinate system yields 

6u = x(~) (3-4) 

whence 



6u being the relative displacement of the end q relative to the end p 
in the corotational system and £ being the strain. 


3.2 FRAME ELEMENT [15] 

A frame element is a structural component which is initially 
straight and which undergoes axial, bending and torsional deformations 
resulting from finite displacements and rotations of its ends. A frame 
element of general cross-section resists loads mainly by three types of 
stresses namely a^, and T xz - For such a frame element, a general 

treatment of shear deformations is quite complex hence only thin-walled 
cross-section frame elements, wherein the plastic strain energy due to 
shear deformations can be ignored, are implemented in ACTION. 

For the linear elastic case a closed form expression for the elas- 
tic strain energy density is developed. 
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For the inelastic case, however, the calculation of the elastic portion 
of the total strain energy density is based on the shear flow theory 
for beams with thin-walled cross-section and the incremental dissipa- 
tive strain energy density U* is assumed to be that due to the effects 
of normal stress and strain alone. 

3.2.1 Geometry of Deformation of a Frame Element 

Figure 3-2 shows the initial position, of a typical frame element, 
specified by the co-ordinates of its end points (the p ane q joints it 
connects) with respect to a fixed system of global axes X, Y and Z. Roll 
orientation is specified by the angle between a reference axis of the 
cross-section and the plane formed by the member and its projection in 
the XY plane. This angle is denoted by y. (If the member is parallel 
to the Z axis, y is measured from the X axis.) 

The longitudinal axis of the element and two reference axes of the 
cross-section form member axes x^, y^ and z^. A vector, {V}, known in 
the global co-ordinate system can be described with respect to the member 
co-ordinate system as vector {V^} given by 

{ V l> = [^{V} (3-6) 

where [T^] is an orthogonal transformation matrix with the property 

[ipVp = [I] (3-7) 

where [I] is the identity matrix. Premultiplication of both sides of 
Eq. (3-6) thus yields 

{V} = [T^ 1 {V 1 > (3-8) 

The transformation matrix, [T^], which describes large angular 
rotations from global axes X, Y, Z to axes x^, y^, z^ are given by 
[16] as 
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Fig. 3-2 Initial Orientation of the Frame Element. 
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[T x ] - 


c c 
y z 


c s 

y z 


"V 


-cs+ssc cc+sss sc 
x z xyz xz xyz xy 


s s + c s c 
■- x z xyz 


-s c +css cc 
x z xyz x y-* 


(3-9) 


where c. = cos<J). and s. = sincj). for i = x, y, z. Angles <f) , $ , (p 
3_ 3_3_ i x y z 

are defined in Fig. 3-2 and the rotations have been performed in the 


order <j) , <f> and <p • From this figure it can be deduced that 
z y x 


cos<P x = cosy 


sinc}^ = siny 


COS(j) = — 

y L 


z - z 

sin<}) ^ — r — 

y L 


cos<i) = 
z 


X - X 

= _a EL 


sin 4> z - 


y - y 
= _3 _R 


where 

l = V(x n ) 2 + (y - y ) 2 and L = \fs/ + (z - z ) 2 

q p q p q P 

Motion of the frame element is characterized as the superposition 

of a rigid body motion and deformations. There is no limit on the size 
of the rigid-body motion; however, the deformations are assumed to be 
small. Motion of the frame element is a function of the three transla- 
tional displacements and the three rotational displacements of the two 
joints which the member connects. 

To define joint displacements and rotations a set of joint axes, at 

the joint p which are initially parallel to the global axes X, Y, Z, 

are introduced. The joint displacements are denoted by the vector IJ 

which has components U, V and W in the X, Y and Z directions respectively. 

The rotations of the joint p, 0 , 0 , 0 , are about the joint axes which 

J r x y z J 

are initially parallel to the global axes. Motion of the element is, 
thus, separated into two parts: a rigid-body motion which is described 

by the displacements and rotations of joint p, and a deformation which is 
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described by the motion of joint q relative to joint p. 

Initially the joint axes are situated parallel to the global axes 
with the origin at joint p. The rigid body motion translates and rotates 
the joint axes to the position x' , y’ , z' as shown in Fig. (3-3). The 
translation is given by the vector and the rotation is described by 
the transformation [T^^ which relates vectors in global system to vec- 
tors in the displaced x' , y’ , z' position of the joint axes, 

{V’> = [T 9 ] {V. } (3-10) 

2 p 1 

The transformation [T„] is identical in form to [T 1 ] of Eq. (3-6) ex- 

2 p 1 

cept that the angles 0 ,0 ,0 are used in place of <J> , <j) , (j) . 

^ & xp yp zp r x’ y T z 

The rotations must be specified in the order 0,0,0 to be able 

r zp yp xp 

to use a form for [ T^ ] similar to [ ] . The member axes remain fixed 
with respect to the joint axes. Rigid-body motion carries the joint 
axes from their initial position parallel to the global axes) to the 
position x' , y', z’ . Likewise, the member axes x^, y^, z^ become axes 
x^y y 2 , z 2 a f fcer the rigid-body motion. However, the new member axes 
X 2 ’ ^2* Z 2 are or: ’'- en ted exactly in the same manner with respect to the 
joint axes x’ , y' , z 1 as were the member axes x^, y^, z^ with respect 
to the global axes. Hence, any vector (V 2 ) described with respect to 
the new member axes x 2 , y 2 > z 2 re l at ed to the vector {V } described 
with respect to the x’ , y ? , z' axes by the relation 


{V 2 > = [T 1 ] {V* } 

Accordingly, use of Eq. (3-9) yields 
{V 2 } = [Tj^] [T 2 ] p {V 1 } 


(3-11) 


{ V = [ T 3 ] p {V l } 


where 


(3-12) 
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(3-13) 


[T 3 ] P = [T 1 ]IT 2 ] p 

Since, both [T^J and [T^^ are orthogonal transformations, their pro- 
duct is also an orthogonal transformation i.e. 

'Vp'Vp - [I i 

Thus, vectors specified in the global axes can be expressed in the de- 


formation axes (member axes after rigid body motion) x 2 » y 2 » z 2 by the 
relation 

{V-} = [T„ ] T {V 0 } (3-14) 

± J p 2 

To describe the deformation of the element, the displacements of 


the joint q need be expressed with respect to the deformation axes. 
From Fig. 3-3 it follows 


6 = R + U - L- U - R 

u -q -q - p -p 


= (R " * > " L + (U - U ). 

— q P — — q — p 

where 6u. is a vector described with respect to the deformation axes. 
Use of Eq. (3-12) then yields 


6u j 

i / X - X i 

i t 

L ' 


ru - u x 


| l q p 

| 1 

j 


q p j 

6v 



° ! 

+ [T„ ] 

v - v ( 


3 p ) q p 

\ J 


P 

q p/ 

6w 

^ ( Z - Z ! 

n n ' 

1 1 

, 0 ] 

1 | 

w - w ) 

n D 7 


(3-15) 


where <5u, 6v, 6w are the displacements of joint q relative to joint p 
along x 2> y 2 , z 2 axes respectively; X, Y, Z are the co-ordinates and U, 
V, W the translational displacements of the joints measured with respect 
to the global axes. 

For numerical calculations, the arrangement of terms in Eq. (3-15) 
is poor, requiring the differencing of nearly equal numbers. Since 
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The terms are rearranged to: 


= [T x ] 


[T,] 

4 p 


AU 

+ [T 0 ]_ <[AV 

AW' 


2 p 


(3-16) 


where [T^ = [T 2 ] - [X] 


(3-17) 

and A is a difference operator for q and p end values in the global 

axes. After manipulation of the trigonometric terms in [T ? ] , [T,] 

AT) Q "P 


may be expressed as 


[ Vp " 


r -2(S^ C 2 +C 2 S 2 ) 
A y2 L z2 y2 b z2 j 


-c s +s s c 

x z x y z 


. S S +C S C 
L_ x z x y z 


C S 

y z 


-S 


-2(S 2 0 C 2 0 +C 2 o S 2 0 )+S S S SC 
x2 z2 x2 z2 x y z x y 


-S C +C S S 
x z x y z 


~ 2 (S x2 C y2 +C x2 S y2^J 


(3-18) 


where the subscript "2" denotes one-half of the indicated angle, the 

angles involved being 0 , 0 and 0 

° xp yp zp 

To complete the description of the deformation, expressions for rota- 
tion of joint q with respect to the deformation axes are required. In 
the development of these expressions, the rotations of each joint are 
permitted to be large but the differences between the rotations defined 
by 


A0=0 -e,A0=e -e,A0=e -e, 

x xq xp y yq yp z zq zp 
are assumed to be small, so that cosA0~l and sinA0-A0. The intent is 

to permit large rigid body motion of the element but limit the defor- 
mation of the element to small rotations. With this restriction, the 
relative rotations \ lj) and \p^ of the joint q with respect to the 
deformation axes are similarly given by 





P 



[TJ 
3 p 



(3-19) 
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Equations (3-19) have been obtained based on the fact that infinitesimal 
rotations can be treated as a vector. 

Thus, deformation of the frame element is specified by the rela- 
tive displacements and rotations of joint q with respect to joint p. 

The relative displacements 6u, 6v, Sw are given by Eqs. (3-16) and the 
small rotations are given by Eqs. (3-19). These equations 

contain large angle transformations which entail considerable calcula- 
tion effort. For problems with small joint displacements and rotations, 
approximate transformations may be 1 employed with considerable savings 
in calculations. This is achieved by replacing the trigonometric func- 
tions by their power series expansion and retaining only terms of the 
second-order. The resulting approximation is given by: 


Sv 

. 6w 


[T^ 


r- ¥ e y + e ? 


- e +ee 

z x y 


- e 


i_ 


L e + 9 0 

y x z 


12 2 

f(e A + ep 

2 x z 


0 +00 
x y z 


12 2 

4(0 + e XJ 

2 x y 


r 1 

-0 


0 - 0 ' 
z y 


L e -0 

y x 


x 

1J 


PI 

(AW) 


(*x) 

/ A0 

CD 

< 

CD 

1 

1 x 

yp z i 

x > - [TJ 

{ A0 

+ 0 A0 ) 

1 

) y 

xp z l 

(-0 

A0 + A0 ' 


xp y 


(3-20a) 


(3-20b 


If all second-order terms are neglected in the above expressions, the 
linear deformation equations are obtained 
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(3-21a) 



(3-21b) 


These results may be used for problems where the displacements and rota- 

2 

tions are small (i.e., (rotations) << order (relative joint displace- 
ments)). If non-linear deformation effects are to be represented as in 
the case of buckling problems, the second-order approximation must be 
employed. 

The parameters 6u, 6v, <5w, ^ , \p z are the generalized displace- 

ments as seen at the end point q relative to the end point p of the ele- 
ment. Deformation along the length of the element is described by intro- 
ducing displacement functions u q (x), v^(x), w^(x) of a longitudinal 
reference axis and the angle of twist 3 q (x) about the reference axis as 
shown in Fig. 3-4. The x, y, z axes of Fig. 3-4 are the > z 2 

axes of Fig. 3-3. The subscript is dropped to simplify the notation in 
the development to follow. 

The displacement functions u q (x), v q (x), w q (x) and $ q (x) define 
deformations of the reference axis. Displacements of points off the 
reference axis are obtained by making the usual kinematic assumptions 
of the engineering theory of beams under bending and torsion. For the 
sake of simplicity, the equations of equilibrium with bending and tor- 
sional deformations, lateral displacements and twists are referenced 
to a longitudinal axis through the shear center and are denoted by v(x) , 
w(x) and B(x) respectively as shown in Fig. 3-4. Based on the results 
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Fig. 3-4 Generalized Displacements of the Frame Element. 


2 } 


from strength of materials and elasticity [17] the following expressions 
are assumed for the displacements u*, v*, w* of points off the centroidal 
reference axis: 


u*(x) = u (x) 


dv, (x) dw, (x) dv (x) dw (x) 

y ~ Z + ^(y.z) — gj— + 4> 2 (y,z) — Ir - 


+ 4> 3 (y»z) 


dg(x) 

dx 


v*(x) = v,(x) + v (x) -(z-z )B(x) 

D s s 

w*(x) = w b (x) + w g (x) + (y-y s ) 3(x> 

(3-22) 

The subscripts "b" and "s" indicate lateral displacement due to bending 
and shear and the functions <J>^, <J> 3 describe cross-section warping. 

The total lateral displacement at the shear center is the sum of the 
bending and shear components 

v = v, + v , w = w, + w (3-23) 

b s b s 

The deformations of the reference axis can be obtained by setting y = 
z = 0 in Eqs. (3-22). It is implicitly assumed in Eqs. (3-22) that lon- 
gitudinal warping is unrestrained and that plane sections remain plane 
during stretching and bending even though the total deformation produces 
nonplanar cross-sections. 

For thin-walled closed cross-sections, the stresses produced by re- 
strained warping are small. The same is not true of a thin-walled open 
section where a deformation model with restrained warping would be de- 
sirable. The assumptions of plane sections remaining plane during bend- 
ing deformations and of plane rotations during torsional deformations 
are quite adequate for inelastic material behavior. However, warping 
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functions become complex functions of the distortion parameters and the 
material properties. The secondary, self-equilibrating stresses due to 
restrained warping, if any, are neglected for elastic and inelastic 
material behavior for sake of convenience. 


3.2.2 Strain-Displacement Relations 

With the assumption of unrestrained warping, the axial strain consis- 


tent with the deformations assumed in Eqs. (3-22) is given by 


. du _ A, 
£ x dx ^ 0 


.2 

d w. 


,dv s 2 


dw s 2 


+ t ef) + c~) 


"dx 


dx 


(3-24) 


dx dx 

Implicit in the Eq. (3-24) is the assumption that squares of rotations are 
negligible in comparison to unity. Equation (3-24) is adequate for des- 
cribing the initiation of buckling but not adequate for describing post 
buckling behavior. A rigorous treatment of shear distortion based on 
the theory of elasticity is quite involved even for the elastic case. 
Hence, to simplify the treatment of shear distortion, especially for 
the inelastic range, an approximate strength of materials approach is 
used. 


3.2.3 Defor m ation Modeshapes: Linear-Elastic Material 

Based on the results from engineering beam theory, deformation mode- 
shapes are assumed and axial and shear strains are derived as functions 
of relative joint displacements. The deformation modeshapes describe 
deformation along the length of the member as caused by relative joint 
displacements Su, 6v, <5w, ^ x > 'Py* *P Z defined in Section 3.2.1. Neglect- 
ing shear deformations and assuming unrestrained warping these results 
are 
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u <5u 
L = n T 


V b ,-2 0 3. 6v b , ,3 2 

— = (3n - 2n ) — jj- + (n - n )i> z 


? b 2 0 3 N 6w b ,3 2v , 

— = (3n - 2n ) — - (n - n W 


(3-25 a-d) 


6 = ml* 


where. 


6v = <5v. + z lb , <5w = 6w, - y ip 

b s T x b ^s T x 


(3-25 e,f) 


ri = x/L and y g , z g are the co-ordinates of the shear center of 
the cross-section of the beam. 

The stress resultants associated with the deformation shapes of Eqs. (3-25) 
are given by 


N = EA 


du 

dx 


d 2 V 

M = El 7T 

z z , 2 

dx 


,2 

d w. 


El 


yz dx 2 


A 2 A 2 

d v d w, 

- M = El El tt 

y yz dx 2 y A J 


dx 


(3-26 a-d) 


T = GJ 4^ 

dx 

where, 

O o 

A = / dA, I = / z dA, I = / y dA, I = / yzdA 
A y A Z A yZ A 

the parameter J is the torsional constant for the cross-section and the 

parameters E and G are the material elastic moduli in tension and shear, 

respectively. Equilibrium of the frame element in Fig. (3-5), gives 
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dM 

V = - El 

y dx z 


,3 

d v. 


dx 


~ - El 
3 yz 


,3 

d w. 


dx~ 


dM 


V = — 7 ^ = - El 
z dx yz 


d 3 v d 3 w 

— El 
3 y 


(3-26e , f ) 


dx 3 dx 

In the development of Eqs. (3-25) , it is assumed that the x-axis or 
the straight line defined by the joints p and q, passes through the 
centroid of the cross-section. Thus, eccentricity between the struc- 
tural joints and the centroidal axis of the element is not represented. 
Note that in Eqs. (3-26) torsion and bending are uncoupled, since twist- 
ing and lateral displacements are referenced to the shear center. 

The effect of shear deformation due to bending is approximated. 

The gross shear distortion of the beam is described by an effective 
shear strain. From the geometry of shear deformation, the shear strain- 
displacement relations are given by 


Y„ 


dv 
s 

dx 


Y r 


dw 
S 

dx 


sxy dx sxz 

where the subscript "s" indicates the shear component of deflection. 

The shear strains produce shear forces based on a gross response de- 
scribed by the expressions 


(3-27a,b) 


dv dw 

_ AG s AG s 

y k dx k dx 

y yz 


(3-28a) 


V 

z 


AG 


k 

zy 


dv 
s^ 

dx 



dw 
s_ 

dx 


(3-28b) 


The constants k , k , k , k are to be determined so that the gross 
y yz zy z 

shear model best represents the shear deformation in bending. From 


among the various methods available for the calculation of the shear 
constants the maximum-strain method is the one used in the ACTION simu- 


lator. In this method, the maximum strain at the centroidal plane is 
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assumed to be the effective strain and 


k = 

y 


k = 
z 


(Y ) 

xy centroid 


(V y /AG) 


Q V 

z y . 


b I G 
z z j 


„ . . , (V /AG) = 

centroid / y 


<Y W > 


xz centroid 
(V z /AG) 


Q V 
y z 


b I G, 

L yy J 


centroid / <V AG > 


PH 
■ PH 


centroid 
(3-29 a,b) 

centroid 


where Q and Q are the first moments of area (on either side of the cen- 
z y 

troidal axis) about the centroidal axes z and y respectively. These 
constants are easily evaluated and results are found in strength of 
materials books [18] for a variety of common cross-sections. 

The constants k y , k yz> k zy , ^ z must satisfy the reciprocity rela- 


tions 


I I 

k = k y - = k = k Z 

yz z I zy y I 

yz J J yz 


(3-29 c, d) 


However, if these constants are obtained by the maximum-strain method 
as outlined above, they may not satisfy the reciprocity relations for 
certain cross-sections. 

Next, solving for v g and W g from Eqs. (3-28) with the help of Eqs. 
(3-26 e,f) and (3-29 c,d) gives 


dv El d 3 v, 

s _ _ , z b 

dx y GA , 3 

dx 

which have the solutions 
6v 


dw 


s s 

L = n — 


s 

dx 


dw 


El dw. 


= - k 


z GA 


dx 


where 


5v 


12k El 

. y . . z 

2 

GAL 


L “ 


6V b 1 , 
IT “2 K 


(3-30 a,b) 


6w 


12k El 

z y 

GAL 2 


<5w. 


I *y] 


(3-30 c , d ) 


Equations (3-30) are simple linear mode shapes that describe the defor- 
mation caused by bending shear. Adding these to the initial mode shapes 
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from Eqs. (3-25) gives 


u 

L “ 11 


6u 

L 


r 9 q 6v H q ? 5v 

J = (3n - 2rr ) -~ + (n - n )ip z + n 


IT + n <“£> 


(3-31 a-d) 


J = (3n 2 - 2n 3 ) - (n 3 - n 2 )i^ y + n 


6w y 

s , s, , 

— - ” ^ K 


6 ‘'K 

2 y s 

where Sv = <Sv, + 6v + (— =r)\Jj , 6w = &w, + 6w - (— )\Jj 

D S Lx D S Lj X 

The deformation modeshapes of Eqs. (3-31) represent a linear defor- 
mation theory. These results may be extended to include the nonlinear 
coupling between axial and lateral deformations. An axial deformation 
of the following form is assumed 

du = K _ I rdv“| 2 + fdw] 2 (3-32a) 

dx K 2 LdxJ LdxJ 

. 

where K is a constant given by the relation 


K = ^ + -T / 

L 2L 0 


ffl + Cf] 


dn 


(3-32b) 


Equation (3-32) replaces the first of Eqs. (3-31). 

Collecting results from Eqs. (3-31), (3-32) the deformation mode 
shapes for the linear elastic frame element are summarized as follows: 


^=^-2l[0 2+( ^ 2 ] 


Sv 


J = (3n 2 - 2n 3 ) -~ + (n 3 - n 2 )^ + n 
l = (3n 2 - 2n 3 ) ^ - (n 3 - n 2 )^ + n ^ 


3 = n $ 


X 


(3-33 a-d) 
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where 


6 = <5v + Sv + z \ p , 6 = 6w, + 6w - y ip ( 3-33 e,f) 

v b s s x w b s J s T x v * ' 

and K is given by Eq. (3-32b) . 

Equations (3-30) and 3-33 e,f) can be solved simultaneously to 
obtain 


6v, 


6v 

L 


a 

-JL 

2 




<5w. 


5w 

L 


a 


, y S , 
+ -r ^ 


1 + a 


1 + a 


where , 


a = 


12k El 

y z 

2 

GAL 


12k El 
z y 

a z = 2 

GAL 


(3-34a,b) 


The parameters and are a measure of the relative importance of 
shear deformation. For a lateral displacement imposed at the end of the 
element, a is the ratio of shear deformation to bending deformation. 

With the deformation mode shapes defined the constant K of Eq. 
(3-32b) can be evaluated. The expression for K is 

r o <5v, o i ^v, i o <5v ^v, i 

K = ^ + lhr> 1 - T5-r* 2 +TJ*z + ~r h r + l^r> 


3 /H,2 ^ 1 H ^ 1 ,2 ^ 5w s 
+ 5 ( L ^ + 10 L ^y + 15 ^y + L 


6 w. 


, . <5w 

_k + I_£) 

L 2 L ' 


(3-35) 


3.2.4 Deformation Mode Shapes: Inelastic Material 

The determination of deformation mode shapes for inelastic action 
in the frame element is not feasible in a general analytical formula- 
tion. For approximate treatment of the inelastic case the elastic 
mode shapes of Eqs. (3-33) are used. For consistency of assumptions, 
however, it can be shown [19] that the linear axial deformation mode- 
shape of Eq. (3-25a) has to be replaced by a quadratic variation of 
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axial displacement u along the length of the element. This is achieved 
by the introduction of an additional node, r, at the center of the frame 
element. Only the axial displacement u^. of this node is monitored. 

Thus 


u 6u- 6u 

r = 371 tt * n 


2 2 ^ U 1 2 ^ u 2 
_2 _ 2n 2 _1 + 2n _ 


(3-36) 


where 

6u- = u - u 

1 r p 

and 5u n = u - u . 

2 q r 

The transformation relation between the global displacement AU^ and 
6u^ can be shown to be 


5u, 


where 


A 1 = T 


A 2 = T 


A 3 = T 


1 

AU 1 

r 1 

- A 

6v* 


A - 

5w* 

A 1 

1 L 


A 2 

L 


A 3 

L 

,11 

T 11 

4_ 

T 21 

T 12 

4- 

T 31 

T 13 

2p 

T 1 

i 

T 2p 

T 1 

T 

T 2 P 

T 1 

,11 

T 21 

4- 

T 21 

T 22 

4- 

T 31 

T 23 

2p 

T 1 

T 

2P 

T 1 

i 

T 2 P 

T 1 

,11 

T 31 

_u 

T 21 

T 32 

_L 

T 31 

t 33 

2p 

T 1 

I 

2p 

T 1 

T 

T 2p 

T 1 


+ 2L (A 4 + A 5 + A 6 )] 


(3-37) 


A 4 = (T ^ Ax + T ip Ay + T 4 X p 3 Az)T 2p 

A 5 " (T 4 P X Ax + T 4p Ay + T 4p Az)T 2p 

A 6 = (T 4p AX + T 4p Ay + T 4p Az)T 2p 

6v* JL r ^ V b . 1 . JL _1 . 

L 2 ( L ) ‘8^z + 2L + 2s^x 


(3-38a-h) 


6w* 1 , 6 V 1 1 6w s 

~ = 2 ( X )+ 8^y + 


2 L 


1 i 

2 y s lJJ x 


and T. 1 ^ is the i-jth element of the matrix T, of Section 3.2.1. Equation 
k k 
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(3-37) assumes that is nonzero. If this is not the case, similar 
equations can be derived in terms of AV^ or AW^. 

It seems appropriate to remark in passing that this feature of the 
frame element using a quadratic variation for the axial displacement 
field is not necessary for the linear case wherein the deformations 
are referenced with respect to the centroidal axis. However, if used 
with deformations being referenced with respect to the centroidal axes, 
the quadratic variation degenerates to a linear variation. More impor- 
tantly, the feature could be exploited to analyze linear response using 
reference axes which are not coincident with the centroidal axes and/or 
for the purposes of simulating rigid links. Of course, this then implies 
that the strain energy of deformations cannot be computed using the usual 
closed form expression of the linear elastic case but rather be computed 
using numerical integration as in the inelastic case to be described 
in Section 6.4. 

For inelastic response, the relative contributions of bending and 
shear to the total lateral deformations are initially unknown. An 
iterative solution for the magnitude of shear deformation is required. The 
iterative solution is implemented by adjusting the relative contribu- 
tions of bending and shear iteratively until the shear deformation as 
measured by the effective shear strain, converges. The following steps 
are involved: 

1. Initial magnitudes of shear and bending deformation are assumed 
based on the elastic solutions of Eqs. (3-30c,d) and (3-34). 

2. An analysis is performed for the stresses and strains in the 

element. The effective shear strains Y , y are evaluated and Eqs. 

' sxy sxz 

(3-27a,b) and (3-30a,b) are used to obtain new measures of the shear 


31 


deformation components 


(<5v ) = y L 

s new sxy 


(<$w ) = y L 

s new sxz 


(3-39a,b) 

3. The new shear deformations are compared to the old, normalized 
by a measure of the total deformations from Eqs. (3-30c,d). 


|(6v ) - (6v ) _ ,| 

1 s new s old 1 

|6v | + |6v I + ]!/ 2 ?TI 


| ( 6w ) - ( Sw ) , 

1 s new s old 

"|6w, [ + l6w| + |l/2 f L| 


« 1 


« 1 


(3-40a,b) 


If these ratios are sufficiently small (equal to 0.10) the iteration is 
stopped, if not, the iteration is repeated from step "2" with the new 
values of shear deformation. 

Consideration of shear deformation requires a considerable cal- 
culation effort for the inelastic case. Therefore, shear deformation 
should be considered only for those elements where it is judged to be 
of significance. 


3.2.5 Shear Flow Theory : 

Figure 3-6 shows the forces acting on an element of a thin-walled 
frame member. With the usual assumptions of shear flow theory of thin- 
walled members [ 18] equilibrium of forces in the longitudinal direction 
yields 


= - $ tds 

o dx 


(3-41) 


where q^ is the shear flow at the origin of s. Equation (3-41) gives the 
variation of shear flow around the perimeter and along the length of 
the beam as a function of the longitudinal stress Q. 

Subsequent treatment of shear stress by the shear flow theory varies 
depending upon whether the cross-section is open or closed. For thin- 
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walled open cross section, the. torsional shear stresses vary while the 
shear stresses due to bending shear are constant across the thickness. 
For thin-walled closed sections the torsional and bending stresses are 
constant across the wall thickness hence both are described by the shear 
flow theory. Thus, in the case of the IE section element (a thin-walled 
open section) the shear stresses due to torsion are neglected both for 
the elastic and inelastic cases, although for the elastic case the 
strain energy due to torsion is accounted for. In this regard the 
treatment of the IE section element is inconsistent. 

For closed cross-sections, Eq. (3-41) describes the change in q 
around the perimeter. The integration constant in Eq. (3-41) is 
selected so that the shear strain, when integrated around the perimeter 
gives the prescribed twist i.e. 


d3 ; 

dx 2A, 


yds 


(3-42) 


where Aq is the area enclosed by the section and y is the shear strain. 

Since, in the formulation of the material model of Section 4, it is 
assumed that the shear response is elastic. In this case Eq. (3-42) 
becomes 


A6 = , 1 , A 3. ds 

dx 2A„G J t Qb 


(3-43) 


0 

Substitution for q from Eq. (3-41) in the above equation yields 

s r q. 


f - / 0 £ c ds]ds + ^ds 


or 


<ln = 


2 A o G 


M _<6il 

dx__J___t 


ds 


(3-44a) 


ds 


where 
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(3-44b) 


0 




It can be shown that the denominator of Eq. (3-44a) is given by 
2 

ds = — — (3-45) 


$ A ■ 4A ° 


where J is the torsional constant for the section. 

If q (s) denotes the total shear flow due to bending and torsion, 

q^(s) the shear flow due to bending alone and q fc (s) the shear flow due 

to torsion then the latter is given by 

q (s) = — ^q ds (3-46) 

t s 0 ^ 

where s^ is the perimeter of the section. Hence 

q fe (s) = q(s) - q fc (s) (3-47) 

3 . 3 MEMBRANE ELEMENT 

A membrane element is a plane triangular thin element which under- 
goes only in-plane deformations resulting from displacements of its 
vertices. No out-of-plane deformations are admitted, but the element 
can undergo large rigid body motions. The usual assumptions of the 
engineering theory of membranes are implied in the development of the 
element stiffness properties. 


3.3.1 Geometry of Deformation of a Membrane Element , 

Figure 3-7 shows the initial and deformed positions of a typical 
membrane element. The initial undeformed position of the triangular 
element is specified by the co-ordinates of its three vertices, p°, 
q ° and r° (taken in a counterclockwise sense) with respect to the fixed 
system of global axes X, Y and Z. The vertices of this triangle after 
deformation are labeled as p, q and r. Local x 1 and y' axes are chosen 
to lie in the plane of the deformed triangle with the x’-axis coinciding 
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Fig. 3-7 Deformation of the Membrane Element. 
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with the side pq. The z f -axis is accordingly perpendicular to the plane 
of the deformed triangle. a and 3 denote the angles at the vertices p° 
and p before and after deformation respectively. If IL, and (i = 
p°, q°, r°) are the displacement components of the ith vertex referenced 
with respect to the X, Y and Z axes respectively, then the co-ordinates 
of the vertices, p, q and r with respect to these axes are 

x. = X. + u.; y. = Y. + V., z. = Z. + w. (i = p,q,r). (3-48) 

i i i u i i xix 


Unit vectors e-^ and e along the sides pq and pr are given by 

AAA 

x i + y j + z k 


2 = -HE 

1 


_re: 


qp = (x i + y j+z k)/R 
/x 2 ' ' + y z ~ + z 7- ‘ qp qp qp 


qp 


qp 


qp 


(3-49 a, c) 


xi + yj+zk /n /s 

- rp - rp ■■■ ■a- - (X. 1 + y 1 + 


/x 2 + y 2 + z z 

rp rp rp 


rp 


rp* 


z k)/Q 
rp 


where 


x.. = x.-x., y.. = y.-y., z. . = z.-z. 
ij i 1 1 3 i J ij iJ 

Accordingly, unit vector e^ normal to the plane of the triangle is 

e, x e [ (y z - y z )i+(z x -z x )j+(x y -x y )k] 

- = = q p r p 3 rp qp qp rp rp qp qp rp rp^qp 

3 sin$ D 

(3-50a) 

where 

D = [ (y z -y z ) 2 +(z x -z x ) 2 +(x y -x y ^2^1/2 (3-50b) 
qp rp 'rp qp qp rp rp qp qp'rp rp qp 

The unit vector e 2 (along the y’-axis) is then given by 


e 2 = e 3 x e;L 


_1 

DR 


{x (y 2 +z 2 ) - x (y y +z z ) }i 
rp qp qp qp qp rp qp r P y 

+{y (z 2 +x 2 ) - y (z z +x x ) } j 

rp v qp qp' 'qp ' qp rp qp rp' 

+{z (x 2 +y 2 ) - z (x x +y y )}k 

rp qp ' qp' qp qp rp 3 qp 3 rp _ 


(3-51) 
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Equations (3-49) through (3-51) together yield the required transforma- 
tion matrix between global and local co-ordinate systems, i.e.. 


or 
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x y z 

\ = _a£ > = _2R > = 

A 1 R * X 2 R * j R 
2 2 

VU = [x (y +z ) - x (yy+zz) /DR 

1 rp qp qp qp qp rp qp rp 

y 0 = [y (z 2 + x 2 ) - y (z z +x x )]/DR 

2 LJ rp qp qp 'qp qp rp qp rp' 

2 2 

= [z (x +y ) - z (x x +y y )]/DR 

3 rp qp 'qp qp v qp rp 'qp'rp /J ' 

V = (y z -y z )/D, V = (z x -z x ) /D 
1 'qp rp 'rp qp 2 qp rp rp qp 


(3-52a) 


(3-52b) 


(3-52c,g) 


and = (x y -x y )/D 
3 qp' rp rp'qp 

3.3.2 Deformation Mode Shapes, Stresses and Strains 

The simplest form of a plane stress triangular element is the con- 
stant strain element. The assumed displacement field as a function of 
the local co-ordinates is 


u = a + a n x + a 0 y 
o 1 2 

v = b + b n x + b„y 
o 1 2 J 


(3-53a,b) 


The constants in Eq. (3-53) are evaluated from the conditions 
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(3-54a,c) 


u(0,0) = u , v(0,0) = v 
P P 

u(R°,0) = , v(R°,0) = 

u(Q°cosa, Q°sina) = u^ , v(Q°cosa, Q°sina) = v 
Substitution of Eqs. (3-54) into (3-53) yields 
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6u 


u = u + ( ~ o ) x + [ 775 — ; 

p v R Q sina 


Sv 
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V = V P + ( 1T )X + [ Q^sina “ 1* Cota] y 


(3-55a,b) 


where 


6u = u - u , 

q q p 

<5v = v - v 

q q P 

(3-55c) 

6u = u - u , 
r r p 

<5v = v - v 

r r p 

(3-55d) 


The rigid body motion of the element in the local co-ordinate system is 

eliminated by setting u = v = v =0. The strains £ , £ and y 

p p q xx yy xy 

for finite deformations are the components of the Green's strain tensor 
which for the two dimensional case reduce to: 
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(3-56a) 


(3-56b) 


(3-56c) 
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From Fig. 3-15 it is evident that 

<5u = (R — R°); <$u = (Qcos3 - Q^cosct) ; <5v = (Qsin3 - Q°sina) (3-57) 

^ r r 

where 


cosa = (X X +YY +ZZ ) /Q°R C 
qp rp qp rp qp rp 7 H 


cos3 =(xx + y y + z z )/QR 
qp rp qp rp qp rp 


(3-58a) 
(3-5 8b) 

Equations (3-56) through (3-58) complete the description of the strain 
field as a function of the relative nodal displacements which in turn 
are functions of the global displacements. It should be noted that non- 
linear terms in the strain displacement relations, Eqs. (3-56) imply 
that the formulation admits arbitrarily large rotations but at best 
moderately large strains. 

3. 3. 2.1 Stresses and Nodal Forces of the Elastic Membrane Element : The 

stresses in the constant strain membrane element, assumed to occur at its 


centroid, are given by 
E 


xx 


(1-V ) 


■5 (e + vg ) ; o 
2 X xx 


= tt - (e + ve ) : 

yy yy (i_ v 2 ) yy xx 


T xy 2 (1+v) Y xy (3-59) 

The corresponding nodal loads in local co-ordinates can then be obtained 


by the relations [20] 
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Section 4 


MATERIAL MODEL 

4.1 Basic Assumptions 

This section provides the theoretical basis along with its under- 
lying assumptions for the material behavior in the inelastic range. 

By material behavior, we imply the prediction of the stress state and 
the strain energy density at a point in a material for a known strain 
state. Two theories are available for describing the material behavior 
in the inelastic range: (i) the deformation or total-strain theory and 

(ii) the flow or incremental strain theory. The difference between the 
two theories lies in the fact that the deformations predicted for the 
volume element by the former theory are independent of the loading path 
while those predicted by the latter theory are path dependent. Further- 
more, the deformation or total strain theory postulates the existence 
of a strain energy density function in terras of total strains while the 
flow or incremental theory postulates the existence of an incremental 
strain energy density function in terms of the incremental strains. In- 
cremental strain theory has the advantage that it describes more fully 
the behavior of a volume element but it has the disadvantage that the 
thoery requires time consuming numerical analysis. The total strain 
theory on the other hand has the advantage of mathematical simplicity 
but does not conform to physical reality for some problems [21]. Be- 
cause of its mathematical simplicity ACTION uses the deformation or total- 
strain theory. It is believed, however, that even at the expense of 
the complexity of the material model there may be a slight computational 
advantage, in terms of the performance of the solution algorithm, to be 
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gained by the use of the incremental strain theory [22]. 

Both the total strain and the incremental flow theories assume that 
the total strain or the strain increment, as the case may be, can be 
decomposed additively into an elastic and a plastic component. Lee [23] 
has shown that such a decomposition is in general not valid for large 
strains but may be justified if the plastic strains are predominant. 
Incidently, Hencky’s total strain theory assumes that in the strain 
hardening range the inelastic component of the total strain is predomin- 
ant. Thus, in view of [23], of the two theories Hencky’s total strain 
theory would appear to be more consistent for cases wherein the strains 
are large. 

Hencky’s theory [21] embodies three hypotheses: (i) the principal 

axes of stress and strain coincide; (ii) Mohr's circle diagram of stress 

and strain are similar at any stage in the inelastic deformation; and 

(iii) volume changes are elastic i.e.> the inelastic deformations are in- 
compressible or that Vp = j-. These hypotheses lead to the following 
relation 

= -^L = constant (4-1) 

°xx-°yy 2 t 

xy 

Furthermore, Hencky's total strain theory assumes that the strain energy 
density, W, is a function of the effective strain, £, defined as 
2 2 2 1 2 1/2 

— (e +e +£ e 4- — y ) for a two dimensional stress state 
jzr xx yy xx yy 2 xy 

e = ^ (4-2) 

e for a uniaxial stress state 
xx 

such that 


dW 

de 


a 


22 2 1/2 

(a +a -a a +3x ) for a two dimensional stress state 

xx yy xx yy xy 

(4-3) 

0 for a uniaxial stress state, 
xx 
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The stress strain relations 


3W 3W , 3W 

w = O , -5 = a and -r = 2 t 

3e xx 3e____ yy 3y xy 


(4-4) 


yy 'xy 

obtained from Eq. (4-2) satisfy Eq. (4-1). 

Equations (4-2) and (4-3) imply the use of Hencky’s total strain 
theory along with its assumption that in the strain hardening range 
the inelastic component of the total strain is predominant [21], This 
in a way is consistent with the assumption that the total strain can be 
decomposed into an elastic and a plastic part especially in cases where 
the strains are large [23]. According to reference [23] it is only when 
the plastic strains are predominant that such a decomposition is justi- 
fied. Thus, the present formulation appears to be consistent in the 
strain-hardening range. 

Next, Von Mises criterion is used to predict yielding. Accord- 
ing to this criterion yielding is assumed to occur when the effective 
stress, a, which incidently is also the second invariant of the stress 
tensor or equivalently the octahedral shear stress, reaches the value 


0 ^, the yield point of a uniaxial tension test. 


4 . 2 Modeling of stress-strain curve and the treatment of stress-strain 
history [15] 

The stress-strain curve of the material under uniaxial tension and 

compression is used as the effective stress-effective strain curve. This 

curve is modeled by eight straight line segments as shown in Fig. 4-1. 

Four of the eight segments describe the tensile side of the curve and 

the remaining four describe the compressive side. The left and right hand 

ends of the curve represent material failure. The sign of the quantity 

(e +£ ) i.e., the sign of the first invariant of the strain tensor is 

xx yy 7 
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used for the selection of the tensile or the compressive branch of the 
curve. It is assumed, however, that the material has the same modulus 
in tension and compression. The model permits situations involving 
initial stress and strain. It is required that all points on the curve 
be uniquely determined by the strain parameter and certain indeterminate 
conditions requiring special treatment are not considered. 

From an unloaded state the loading for cyclic increasing or de- 
creasing strain follows the initial linear elastic portion of the curve 
as long as the effective stress does not exceed the stresses correspond- 
ing to tensile and compressive yield points. For increasing effective 
strain beyond the yield points, the loading follows the effective stress- 
strain curve with yielding occurring as described in the previous sec- 
tion. Once into the plastic region, when the effective strain starts 
to decrease, unloading occurs. Unloading in the plastic region is as- 
sumed to be elastic with the load path described by the shape of the 
initial elastic portion of the curve. Points f b' and ' d' in Fig. 4-2 
denote the new tensile and compressive yield points. Loading for addi- 
tional decreasing and increasing strain follows the new elastic curve as 
long as the effective stress does not once again exceed the stresses corres- 
ponding to the new yield points. For straining beyond the new yield 
points, the subsequent loading and unloading is treated as described 
earlier except that the portion of the stress-strain curve between points 
'^q 1 , f d ' in Fig. 4-2 is for all practical purposes forever lost from 
the "memory" of the material. This type of stress-strain behavior gives 
a larger tensile yield stress and a smaller compressive yield stress as a 
result of plastic deformation in tension. The converse is true for plas- 
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Fig. 4-2 Typical Stress-Strain Curve. 
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tic deformation in compression. These changes to the stress-strain 
curve are the strain-hardening and Bauschinger effects. 

Recall that the transient response is based on a stepwise integra- 
tion in time. Each time step corresponds to a loading or an unloading 
increment on the stress-strain curve. Let Gq, correspond to a previous 
time point and G^, represent the state at the current time point. If 
is greater that it is assumed that the material "loads up" only. 

Likewise if is less than Eq, it is assumed that the material loads 
down only. This is illustrated in Fig. 4-3. It must be recognized that 
for any given loading two or more solutions all of which satisfy equil- 
ibrium are possible but all these solutions correspond to different 
strain histories. Hence, for greater than £q the only possible 
solution is the one corresponding to the point T a'. The solution corres- 
ponding to the point T c’ satisfies the constitutive relation but vio- 
lates equilibrium for a member loading up. The point 1 c r is a possible 
solution for greater than £q if and only if the member is first 
loaded up to the point 'b* and loaded down from this point to the point 
’c’. Such two step loading processes are not admitted in the formula- 
tion. 


4.3 Evaluation of Dissipative Strain Energy Density 

For the elastic-plastic response the strain energy density is de- 
composed into an elastic part and an incremental dissipative part. 

Thus, if ’a^’ in Fig. 4-4 denotes the state at some previous time t^ and 
if ? a' denotes the state at time (t^+At) then the incremental dissipa- 
tive strain energy density AW^ and the elastic strain energy density 
are, by assumptions, the appropriate areas under the idealized effective 
stress— effective strain curve as shown in Fig. 4-4. 
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Fig. 4-4 Evaluation of Strain Energy Densities. 


For the truss element is the only strain prescribed and the 
corresponding stress cf^ is assumed constant over the length of the 
element. For the frame element, however, stresses vary throughout the 
volume of the element. Certain number of Lobatto quadrature or refer- 
ence points are used to describe the stress distribution throughout the 
element. At each of these points, the material model is used to deter- 
mine the stress and the strain energy density (W e +AW^) . In the interest 
of simplicity, it is assumed that the shear stress-strain behavior is 
adequately described by linear elastic response. This implies that the 
plastic action is adequately described by the effects of normal stresses 
and strains alone. This approximation should be appropriate for ele- 
ments wherein plastic axial or bending actions predominate. In the case 
of a membrane element all the three stress components are by assumption, 
constant throughout the element and hence as in the case of the truss 
element the stresses and the strain energy density need be evaluated at 
a single point. 

The above ingredients comprise the treatment of loading and un- 
loading in the plastic regions with strain-hardening and Bauschinger ef- 
fects. Although, the model used here is typical of most material be- 
havior it is relatively simple and not all material behavior will fit 
this model. Caution should particularly be exercised when this 
simplified material model is used in conjunction with membrane elements 
(two-dimensional stress states) undergoing cyclic loading. 
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Section 5 


KINEMATIC CONSTRAINTS 

Practical modeling considerations demand that equations of motion 
be augmented by constraint equations. This section describes constraint 
models in ACTION which facilitate representing a rigid tie between joints 
and an impenetrable contact plane. 

5. 1 Rigid Link Element 

A "rigid link" is a finite element which does not deform appre- 
ciably. As the nodes in the structure undergo displacements, the rigid 
element merely translates and rotates without any appreciable deforma- 
tions. Such a link can be used as a connector between two members meet- 
ing at a joint whose ends do not coincide thereby simulating eccentri- 
cally connected members. Rigid links may also be employed advantageously 
when entire sections of the structure undergo very little deformation. 
During these periods the sections can be treated as assemblages of rigid 
links to greatly reduce the number of unknown displacements and thereby 
reduce the computing time. The nodes of a rigid link are located with 
reference to a global co-ordinate system and they are identified as be- 
longing to the same rigid link by input specification to ACTION. Points 
in each rigid link are referenced by a local co-ordinate system originat- 
ing at the primary node of the link. The local co-ordinate system trans- 
lates and rotates with the primary node as deformation takes place. 

Since there is no relative displacement due to deformation in the 
rigid link, the right hand side of Equation (3-15) set equal to zero will 
provide the relationships necessary to describe the motion of the q-end or 
secondary node of the link in terms of the translation components and ro- 
tation components of the p-end or primary node of the link. 


51 



Solving this equation for the global displacement vector of the 
secondary node of the link yields 

{U} = {U} + [T.] T {L} + {R} - {R} (5-1) 

q p 3 J p p q 

This equation is envoked whenever the displacement components of the nodes 
are updated. 

The equation is implemented by searching the list of elements for those 
designated as rigid links. When one is found, the displacement components 
of the secondary nodes are replaced by those calculated by Equation (5-1) . 
Furthermore, given the translation and rotation components of the velocities 
and accelerations of the primary node, the motion of a secondary node is 
determined using simple rigid body kinematics. Finally, knowing the 
kinematics of the secondary nodes the contribution to the total poten- 
tial energy, of inertia forces and loads applied directly to such secon- 
dary nodes, can be determined. 

5 . 2 Impenetrable Contact Plane (Terrain Model) 

The ACTION code includes a model of an impenetrable terrain to simu- 
late a ground plane. The ground plane is assumed to be rigid and flat 
i.e., like a concrete runway. Penetration of the ground is not allowed. 
Resistance to forward motion along the ground plane is provided by 
Coulomb friction. Impact with the plane is represented as a plastic col- 
lision. 

The paragraphs that follow explain the ground plane model in ACTION. 
5.2.1 Node Capture and Release 

The model detects when the nodes of the aircraft model contact the 
ground plane. Once contact has been made these nodes are constrained 
to remain on the ground plane until they are pulled off by the release 
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of internal energy stored within the aircraft model. Thus a node of the 
aircraft model does not bounce off the ground plane due to the effect of 
the coefficient of restitution, but is pulled off when tension would other- 
wise be implied between the structure and the impenetrable plane. 

Nodal displacements are obtained by solving the equations of motion 
for a particular time step. The current position of each node is then 
determined relative to the ground plane. If a node penetrates the ground 
plane, the time step is reduced and integration of the equations of 
motion is repeated. When a node penetrates the ground plane using the 
minimum time step, a check is made to see if its previous position above 
the ground plane was closer to the ground plane than its current posi- 
tion which is beneath the ground plane. (The likelihood of a node fall- 
ing exactly on the ground plane at the end of a time step is remote.) 

If the node’s previous position was closer to the ground plane, integra- 
tion is backed up to the beginning of the current time step and a flag 
is set to signify contact between the node and the ground plane. If 
the node's position at the end of the current time step is closer to the 
ground plane, this position is used to mark the ground plane and a similar 
flag is set. Thus a node is considered to lie on the ground plane when 
it is as close to the ground plane as permitted by the minimum integration 
time step. (In terms of the logic of the code, the ground plane moves 
up or down to coincide with the position of the node.) 

Once contact has been made between a node and the ground plane, it 
must be determined whether or not the node is to be constrained to re- 
main on the ground plane during the next integration time step. This re- 
quires calculation of the total vertical resultant force acting on the 
structural model of the aircraft at the captured node. This force, which 
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is normal to the ground plane and positive upward, is the derivative of 
the strain energy function with respect to a displacement of the node 
in the direction of the outward normal to the ground plane. A positive 
normal force means the structural model is in compression in the vicinty 
of the captured node. In this case the vertical degree of freedom of the 
node is constrained for the next time step, but the node is still free 
to move horizontally along the ground plane. The vertical degree of 
freedom is constrained in the sense that the vertical displacement of the 
node corresponding to its current position on the ground plane is fixed 
for the next time step. A negative normal force means the structural 
model is in tension at the captured node. Then the vertical velocity, 
acceleration and all the higher order time derivatives of displacement 
of the node implied in the temporal algorithm are set equal to zero. 

The free equations of motion provide the basis for predicting subsequent 
behavior of the node until it is recaptured. 

Due to the discretization of time it is assumed that the best approxi- 
mation to the point in time when a negative normal force is found is based 
on the minimum integration time step. Thus if a negative normal force is 
discovered at the end of a non-minimum time step, the time step is reduced 
and the integration is repeated. 

A coulomb friction force is applied to each node captured by the 
plane. The force is applied in the direction opposite to that of the hori- 
zontal velocity components of the captured node. The magnitude of the 
force is equal to the coefficient of friction input by the user times 
the positive resultant nodal force normal to the ground plane. 
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Section 6 


ANALYSIS BY ENERGY MINIMIZATION 
6. 1 Background Information 

It was indicated in Section 2 that two distinct solution approaches 
exist: (i) the vector approach and (ii) the scalar approach. In the 

former approach, the mathematical model is derived on the basis of the 
principle of virtual work and reduced to a system of non-linear second- 
order differential equations in time. In the latter approach, a scalar 
or a potential function associated with the energy of the model is 
introduced, minimization of which yields the desired equilibrium con- 
figuration. In both approaches a temporal finite difference scheme is 
utilized to effectively eliminate time as a variable. As a result, in 
the vector approach the equations of motion are reduced to a system of 
nonlinear algebraic equations in the unknown nodal parameters of the 
finite element model. In the scalar approach which is of relevance to 
this report the problem is reduced to a well known problem in mathema- 
tical programming namely the unconstrained minimization of a nonlinear 
function of several variables. 

For all structural problems with geometric and material nonlinear- 
ities of the type considered herein the required potential function al- 
ways exists. Although, this technique has been hitherto used for mainly 
positive or negative definite systems, other systems which fail to be 
positive or negative definite can be handled by using the least squares 
method or the modified conjugate gradient method with preconditioning 
[24]. In some cases, for such systems displacement incrementation 
rather than load incrementation in conjunction with conventional uncon- 
strained minimization techniques can also be effective [25]. 
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6.2 Solution Basis 


The minimization scheme as applied to the solution of transient 
nonlinear structural analysis problems consists of minimizing a poten- 
tial function associated with the system for an assumed relationship 
between displacements and time. The displacement-time relation for each 
generalized nodal displacement of a finite element model is assumed to 
be of the form [26] 

q e . = B(it) 2 q e . + - B)(At) 2 q 0 . + (At)i 0 . + q 01 (6-la) 


q ei = + C 1 -Y)(At)q 0i + q Qi (6-lb) 

where q ^ is the i-th generalized nodal displacement at the end of the 

• it 

time step and 3 and y are constants. The quantities and can 
now be expressed in terms of the i-th generalized nodal displacement, 
^Oi’ velocity* q 0± and acceleration, q^ at the beginning of the time 
step and the generalized nodal displacement, q £i > at the end of the time 
step. Thus, 
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q ei = q 0i + 


3(At) 


2 ^ q ei “ q 0i' ) “ (At)q 0i " 2 q 0i^ 


(6-2b) 


The equations of equilibrium 
3TT 

M.q . - F. + =0, i = 1,2,. . . ,N (6-3) 

x ex x dq , 
ex 

for an N degree-of-f reedom system with lumped masses can then be written 


M i^0i + 
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2 {(q el - W 
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ex 
(6-3a) 


It can be easily verified that Eqs. (6- 3a) correspond to the necessary 
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conditions for the functional 


S = l ([- 


i=l 


- F. 


23 (At) 


qL - <- 


2 M ei 


3(At) 


2 q 0i + 0(At) q 0i + ( 23 1)q Oi )q ei ]M i 


il(t + At) q ex J+ U + C 


(6-4) 


to be stationary. In Eq. (6-4), U is the strain energy and C is an arbi- 
trary constant. Thus, knowing q^, q Qi and q^ at time t for any given 
load F^ at time (t+At), the functional S may be minimized with respect to 
the generalized nodal displacements, q ^ (i=l,...,N), in order to deter- 
mine the corresponding stable equilibrium configuration. Although, the 
coefficient of in Eq. (6-4) is quadratic in q^, the strain energy, 

U, will in general be a nonlinear function (at the very least a quad- 
ratic for a linear problem) of the generalized nodal displacements q 
(see Section 3). Thus, in general the functional S is highly nonlinear. 
Since, U is a positive semi-definite function for most structural ma- 
terials, the functional S can be seen to be convex. 


6 . 3 Minimization Algorithms 

Of all the algorithms for unconstrained minimization only the quasi- 
Newton or the variable metric algorithms have been more frequently used 
for such nonlinear problems, because of their higher effectiveness [27]. 
Again, unless there exists an algorithm which exploits and maintains 
sparsity of the Hessian approximation during its passage to the minimum, 
one has to resort to some form of a first order conjugate gradient al- 
gorithm for problems wherein N is an extremely large number. 

Beginning with an arbitrary initial guess, these algorithms seek 
a direction of travel and the amount of travel in that direction. The 
manner in which these are sought depends upon the sophistication of 
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the particular algorithm invoked. Most often the directions of travel 
are sought in a manner which guarantees not only a decrease in the value 
of the function but also convergence to the minimum in a finite number 
of iterations (usually N+l for an N dimensional space) in the case of 
quadratic functionals. It is important to note that all functionals in 
question are very nearly quadratic in the neighborhood of the minimum. 
6.3.1 BFGS Variable Metric Algorithm 

The present formulation uses the well-known BFGS (Broyden-Fletcher- 
Goldf arb-Shanno) variable metric algorithm [28] which is supposedly the 
best current variable metric update formula for use in unconstrained 
minimization. This algorithm dispenses with the exact line searches 
while using an update formula which, in the case of a quadratic func- 
tional, guarantees a monotonic convergence of the eigenvalues of the ap- 
proximating matrix to the inverse Hessian. The iterative scheme which 
is begun with the null vector as the initial guess is defined by 

{q} (k+1) - {,}« - a< k W k >{g} (k > (6-5.) 

where 

(g)^ = VS(q^) = gradient of S at q^ k ^ 

(k) 

[H] is a matrix which is designed to approximate in some sense 
the inverse Hessian matrix of S at {q}^ and i s an appropriately 


chosen scalar. The BFGS update formula is given by 

[H] (k+1) - ([I] - A [P])[H] Ck) ([I] - A [P]) + {o}{a} T / g (6-5b) 

where 

[P] = (aHy} T (6-5c) 

3 = {a} T { y} (6-5d) 

(a) = (q} (k+1) - {q} (k) (6-5e) 

(y> = {g} (k+1) - {g} (k) (6-5f ) 
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00 

The a in Eq. (6-5a) is chosen either by a linear search or by a step 
length but so as to maintain the positive definiteness of the matrix 
[H] . The details of the step length calculation may be found in refer- 
ence [29] or reference [30], 

The iterative scheme is begun with an initial guess which for the 
first time step is usually the null vector for {q}, the unknown general- 
ized nodal displacements and the identity matrix for [H], the approxi- 
mation to the inverse hessian. From thereon these quantities at the end 
of the previous time or load step are used as initial guesses for the 
next step and it is this reasonably good approximation to the inverse 
Hessian that is instrumental in giving the second order methods a signi- 
ficant advantage over the first order methods like the conjugate grad- 
ient techniques. The required gradient of S is evaluated analytically 
although doing so by a finite differencing scheme is also possible. The 
use of an analytic gradient however, results in a substantial saving in 
the computational effort. This saving is the result of not only a 
cheaper gradient evaluation but most often a faster convergence of the 
solution because of higher accuracy of all the computed quantities [27], 
The i-th component of the gradient of S can be written as 


3S „ „ . 3U 

t: — — = M.q . - F . 4- -x 

dq . l ei l dq^ 


( 6 - 6 ) 


ei ei 

The term in Eq. (6-6) requiring significant computational effort is 

as it embraces the geometric and material nonlinearities. Thus, 

3q 


ei 


- f f dv = £ / ("wSL.)*, 

k 3q ei k k=l V k di k 3q ei k k 


(6-7) 


M ei k=l 

where 

W = strain energy density; £ and a the effective strain and stress 
effectively as defined in Eqs. (4-2) and (4-3). 
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6.3.2 Powell* s Conjugate Gradient Algorithm 


The algorithm [33] is designed to improve the linear convergence 
rate of the Fletcher- Reeves 1 conjugate gradient algorithm [34] by using 
a restart whose frequency is not the usual n or (n+1) iterations but 
rather one which is dependent on the objective function. Thus the 
basic algorithm may be stated as follows. 

Given {q}, the initial direction of travel {d}, is defined to be 


the steepest descent direction -{g}^ = -{VS}^. For k > 2 

{q>k +1 - {q} k + \{d} k (6-8 a) 

{d> k - -{g} k + S k {d> k _ 1 + Y k {d> t (6-8b) 

6 k " '‘ g i k ((s) k - ^k-l^ d \-l^ s \. " ^k-l^ (6-8c) 

T k = ^ g ^k^ 8 ^t+l * (6-8d) 

where t is initially set equal to one and for k > 2 t is set equal to 
k-1 if 

|{ g }£_l{ g } k l > °.2| |{ g > k l ! 2 - (6-8e) 


in Eq. (6-8a) is determined by a line search which yields not only 
a decrease in the directional derivative but also requires the angle 
between arK ^ to be close to ninety degrees. Next, if the 

inequalities 

-1-2 | |{g} k | | 2 < {d> k {g> k < -0.8| |{g) k | | 2 (6-8f ) 

are not satisfied then {dj-^ is assumed to be not sufficiently downhill 
and the procedure is restarted by letting t = k-1 and redefining {d}^ 
by letting y^ in Eq. (6-8a) to be zero. The procedure is also 
restarted if (k-t) > n by setting t = k-1 and assuming y^ in Eq. (6-8 a) 
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to be zero when k = t+1. For additional details of this algorithm along 
with its fortran listing the interested reader should refer to reference 
[35] which also provides a listing of the subroutine for the BFGS var- 
iable metric algorithm of the previous section. 

The storage requirements of this algorithm exceeds those of 
Fletcher-Reeves ' only slightly but they are certainly very small by 
comparison with those of the variable metric algorithms. Thus, this 
algorithm is intended for extremely large scale problems of the 
order of thousand degrees of freedom or above. However, the 
performance of this algorithm for the solution of such large scale 
problems of relevance to this report remains to be investigated. 

6 . 4 Evaluation of the Function S and its Gradient 

The function S and its gradient must be expressed explicitly or im- 
plicitly as a function of the global generalized nodal displacements of 
the finite element model. 

6.4.1 Function Evaluation 

From a known vector of the generalized nodal variables in the global 
co-ordinate system, consistent with the prescribed boundary conditions, a 
vector of local generalized variables in the co-rotational co-ordinate 
system of each element is established through transformations which are 
functions of its geometry and its rigid body rotations (see Section 3) . 
The assumption of deformation patterns of the element as functions of 
these local generalized nodal variables (interpolating polynomials) 
yields element strains . Recourse to the element material model then 
yields the corresponding stresses and strain energy densities at various 
predetermined quadrature points over the extent of the element. Barring 
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purely elastic response, a simple weighted summation of these quantities 
over the element volume yields stress resultants and strain energies re- 
spectively. For purely elastic response these are provided by well-known 
closed form expressions which can be generated as follows based on the 
results of Section 3. 

6.4. 1.1 Strain Energy Evaluation for Elastic Response 
The total strain energy U may be expressed as 

u - I U (6-9) 

k=l K 

where for purely elastic response a closed form expression for (J^, the 
strain energy of the k-th element, can be obtained from the usual defini- 
tion of the strain energy and the use of the results of Section 3. Out- 
lined below are such expressions for the truss, the frame and the mem- 
brane elements. 


6. 4. 1.1.1 Truss Element 

The strain energy of the k-th truss element is given by 


tJ 


a e dv = 0 
v x x 2 



which from Eq. (3-5) can be written as 

,, _ EAL .DL. 2 

k 2 C 1/ 

wherein E, A and L pertain to the k-th stringer element. 


( 6 - 10 ) 


6. 4. 1.1. 2 Frame Element 

The total strain energy of deformations of finite element due to 
bending, tension and shear is given by 


U k = 


u, , + u, . + u 


kb 


kt 


ks 


where 


( 6 - 11 ) 


62 



U lrU ° l L £ ! dV 


kb 2 J v x 

which from Eqs. (3-24), (3-25) and (3-32) reduces to 

- ^ [r2a + rf { V ( tt >2 - hr>*. + ^l + 2I y Z [( T k)( T k) 

L 

i 6v 6w, . 6w, 0 6w, 

+ 2 <TT *y - IT V - I V. 1 + V O + ( T>y 




u kt = I/ 0 T i d * 


which from Eq. (3-33) reduces to 

/ / GJ . 2 
U kt = 2T 


and 


1 

ti, = ? / (V y + V y )dx 

j^g J J TT 1 CV17 *7 ' C?V7 


0 


y sxy 


which from Eqs. (3-27), (3-28), (3-29) reduces to 
U. 


ks 


r-AT i 6v o i , 6v 6w 

— r— c — ) 2 + i o ---- + 1 ' ' S 

2 L k ^ L ; + yz l k I 

y y z z y 


( 6 - 12 ) 


(6-13a) 


(6-13b) 


6w 


S\ 2- 


k i <-r> + ir ( L > 1 (6 - 13c) 


Frame elements with only five typical cross-sections namely BOX, IE, 
TUBE, ELIP and SORE (Solid Rectangle) are permitted. Explicit express- 
ions for the various section constants associated with each of these 
cross-sections can be found in Appendix A. 


6. 4. 1.4. 3 Membrane Element 

The strain energy of the k-th membrane element can be obtained by 
integrating the strain energy density function W* over the volume of the 
element. The assumption of constant strain within the element simplifies 
this integration to yield 

> M_\0 9 

(6-14a) 


1 1 EAh r 2 . 2 . „ , (1-v) 2 , 

u, = — [e_ + e„__ + 2ve e + y y ] 


-M-, 2 ** yy 


(l-v) 


xx yy 


xy 
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(6-14b) 


where from Figure (3-7) 

A = y Q°R° sin a 

and the strain terms are given by Eqs. (3-56) through (3-58). 


6.4.1. 2 Strain Energy Evaluation for Inelastic Response 

Although closed form analytical expressions for U can be developed 
when the material is elastic the same is not true when the material 
yields. Then the response depends upon the current values of stress 
components and past history. As shown in Section 4 Von Mises 1 yield 
criterion together with Henckey's total strain theory provides a simple 
means of calculating strain energy density distributions throughout an 
element that has yielded. Because total stresses and total strains are 
no longer linearly related recourse must be made to numerical integra- 
tion of the strain energy density over the volume of the element. 

The strain energy density may be decomposed into an elastic part 
and an incremental dissipative part thereby providing an estimate of the 
total energy of the system that has been dissipated through inelastic 
deformations. Thus for a system with m elements 

m . m 

U = strain energy = £ U 1 = £ (J 1 + A(J^ 

i=l i=l G d 


m 


= V ( f W 1 dv + / AWjdv) 

" J IT a •> y (J 

i 


. , v. e 
i=l i 


(6-15) 


where W 1 and AWj are obtained from the material model as described in 
e d 

Section 4 (see Figure 4.3), 

For open or closed cross-sections of the type shown in Fig. (6-1) . 
The integrals in Eq. (6-15) are expressed as sums of integrals over a 
finite number of strips, n^. In keeping with the assumptions of the thin 
walled theory the integrands are assumed constant over the thickness 
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of these strips. Both integrals in Eq. (6-15) are of the form 


n 

s 

l 

j=l 


tj !a 3 


f.dA 

3 


(6-16) 


which using a three point Lobatto quadrature rule in each of the two 
co-ordinate directions in the plane of the strip reduces to 
3 3 n s 


l l l 

1=1 k=l j=l 


t .A.H 0 H. 

3 3 & k 


f j£k 


(6-17) 


and are the weights and f ^ is the value of f at the £-kth 
quadrature point [31] . The exact locations of these quadrature points 
(or equivalently also the stress reference points) for the four dif- 
ferent cross-sections can be found by referring to Appendix A. 


6.4.2 Gradient Evaluation 

The gradient of S as defined by Eq. (6-7) involves the accelera- 
tion vector q . which for a given vector q . , is provided by Eq . (6-12) 


3(At) 


[(q. 


- q 


Oi 


,) - q 0i (At) - (^ - S)q n ,(At) 2 ] 


Oi 


(6-18) 


and the gradient of U with respect to q is given by the general ex- 
pressions in Eqs. (6-8). The expressions in Eqs. (6-8) can, however, 
be simplified and reduced to a closed form if the response is purely 
elastic . 


6. 4. 2.1 Strain Energy Gradient Evaluation for Elastic Response 

For purely elastic response the i-th component of the gradient of 


strain energy can be calculated as 


m r 3li 3r. 

k, , J_- 


3q 


u// ® 

3U - I *r-- I l Ofer 1 ) 


. L. 3q . L m L 3r , 9q 

ex k=l ei k=l j=l j ex 


(6-19) 


where r^ are the local relative degrees of freedom of the q-th node 
relative to the p~th node and n^_ the total of such local relative degrees 
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of freedom. is the strain energy of the k-th member, k =» l,2,...,m 

being the number of elements which have the degree of freedom in 

common. 


6.4. 2. 1.1 Truss Element 


The relative degree of freedom in this case, £ = DL/L, is a func- 
tion of the relative global displacements AU, AV and AW. Thus, 


w- - l < EAL > e ^ 

q x=l q 


From Eqs. (3-2) and (3-3) it follows that 


3e _ (Ax+Au) 
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(6-20) 


(6-21a) 


(6-21b) 


(6-21c) 


6. 4. 2. 1.2 Frame Element [37] 

The relative degrees of freedom in this case are 6u, 6v, , <Sw, , 

b b 

ipy, and ip z and whereas Eqs. (6-12) through (6-14) provide 11 as a 
function of these relative generalized displacements, Eqs. (3-16) through 
(3-21) provide the necessary relationships between them and the relative 


global generalized displacements AU, AV, AW, AG^, A0^ and AG^. 

Thus for the k-th element 

3U, X 

k _ ,ou> 
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(6-22a) 


3 U, 


3(6^) " fl ~ + “y. 1 ^ + K ]> 


<5v. 


<5w. 


12k El . 5v I 
- z - _yz , I_ 1 

2 l k 1 L J 2 ( k I + k I 
L y y z z y 


6w 

X-iT» 


(6-22b) 
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+ — ^ { ( ; J5, + 5-i-) (-J) + 3T ^ } (6 - 22c) 

L 


y z z y 


9U k = GJ 

3ib L V x 
r x 


(6-22d) 


9U 6w 0 1 6v, n 

■ Jf { y O + ! v + 2 Vi ( ^ } - i V 5 

6K El I -a . 6v , <5w 

+ {(-p-xrf + rfx-T^) ♦ ^ <-^» 

y z z y z 


(6-22e) 


9U 


Sv. 


<5w„ 


k - ^ {I r-c-r 1 ) + 4 1 ( 1 .] - 21 [i (-r 1 ) + 4 , U> 


9ii) t 2 lj 'z 

z L 


3 Y z 


yz l 2 L 3 y 


6k El I , . 6w -a <Sv 

-?-* «-f> <ft- + rf > «-r> + r hr” 

y z z y y 


(6-22f ) 


Furthermore, it can be easily verified from Eqs. (3-30c,d), (3-34 a, b) 
that 
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r 96v 
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It remains to derive expressions for the derivatives of the relative 
generalized displacements. For the case of large of deformations pro- 
vided by Eqs. (3-16), (3-17) and (3-19) 

= _ ex ] T [t ] T 

9{U} 1 2 J p U 1 J 

P 
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where the components of [T^J 



are 


T T 21 “ (s x S z + C X S y C Z )(AX + AU) + (- s x c z 

+ c x s y s z > (AY + AV) + c x c y (AZ + AW) 

T T 31 = (c x s z - S x S y C Z )(AX + AU) - (c x C z 
+ s x s y s z > (AY + AV) - s x s y (AZ + AW) 

12 

t t “ - s y c z (AX + AU) - s y s z (AY + AV) 

- c y (AZ + AW) 

22 

= s c c (AX + AU) + s c s (AY + AV) 
i xyz xyz 

- s x Sy (AZ + AW) 

t t 32 = c x c y c z (AX + AU) + c x Cy s z (AY + AV) 

- c x Sy (AZ + AW) 


(6-24c,h) 


(6-25a,h) 
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T y z y z 
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(6-26i) 
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6. 4. 2. 1.3 Membrane Element 

The relative degrees of freedom in this case are $ u r and 

(see Fig. 3-7). Thus for moderately large strains 
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with similar relations for derivatives of x . y and z with respect 

rp rp rp 

to q . . 
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6.4. 2. 2 Strain Energy Gradient Evaluation for Inelastic Response 

For inelastic response the i-th component of the gradient of strain 
energy is evaluated as 
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ei 


where W and £ are defined as in Eqs. (4-4) and (4-5) and the remaining 
equations are defined in Section 6.4. 1.2. As in the case of the strain 
energy, the gradient evaluation for inelastic response necessitates a 
similar integration using quadratures. Except for the change in the 
integrand this evaluation is identical to that described previously. 
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6. 4. 2. 2.1 Truss Element 

In this case £^ = £^ which is also the relative degree of freedom 
(j=l). Thus 

'‘Br . k L 
3 

Expressions for (9r./9q .) are provided by Eqs. (6-21). 
j ei 


6.4. 2. 2. 2 Frame Element 

In this case again £^ = (£ ) . With the aid of Eqs. (3-24) and 


(3-32) this is seen to be 
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where K is given by Eq. (3-35). Thus, for the case of large elastic 
deformations 
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Expressions for (3r../3q ei ) are provided by Eqs. (6-26), (6-27) and 
(6-28). 


When a frame element responds inelastically , an additional node 
is added at the center of the element. The axial displacement field 
is then defined with the aid of two relative axial degrees of freedom 
<5u^ and as evident from Eq. (3-36). Equation (3-37) provides the 
definition of 6u^ in terms of other usual quantities. The derivatives 
of 6u^ and 6u£ can thus be obtained by differentiation of relations 
in Eq. (3-38a-h) . Explicit expressions for such derivatives are however 
too lengthy and cumbersome and have been omitted accordingly. 


6. 4. 2. 2. 3 Membrane Element 

The effective strain, E^ i- n this case is defined by Eq. (4-4). 
Hence, 
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Expressions for (3r./3q .) are provided by Eqs. (6-29). 
J ei 
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Section 7 


SOLUTION ERROR CONTROL 

Errors in the present solution process arise from two sources: 

(i) truncation errors which occur due to truncating the series 
representation of the response Eqs. (6-1) and (ii) round-off errors 
which occur due to the use of a computer with finite digit arithmetic 
precision. In addition to the accuracy considerations of the solution 
process due regard must be taken of its stability. For linear systems 
it has been shown by Goudreau and Taylor [32] that in order to maintain 
unconditional stability of the numerical integration scheme given by 
Eqs. (6-1) the parameter 3 should be chosen such that 

3 > 0.25 (0.5 + y) 2 

with y > 0.50. No such criteria can be postulated for nonlinear 
problems however, with the optimum values of 3 and y being very much 
problem dependent [25]. 

Proper choice of a step size (time or load) is thus crucial. If 
two large a time step is used truncation errors occur due to 
truncating the series representation of the response. Furthermore, 
there is a danger of instability of the numerical integration process. 
If too small a step size is used, computer accuracy limitations arise 
due to limited number of digits used by the computer to represent the 
response. In practice, irrespective of whether a problem is linear 
or not the determination of the optimum size is difficult because of 
the difficulties of accurately measuring the errors due to polynomial 
truncation and computer arithmetic truncation. 

Automatic error control in ACTION consists of varying time steps 
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to assure that the truncation error is tolerable. The errors are 
controlled by adopting time steps so that the equations of motion are 
satisfied to user specified accuracies at midstep times. In other words, 
some norm of the gradient of S in Eq. (6-4) is required to be less than 
a user specified limit. Using Eqs. (6-1) responses are interpolated 
for mid time. An equilibrium check is made using these displacements 
and midtime forces and interpolated accelerations. Similar equilibrium 
checks are also made at the end of the time step. If the errors are 
excessive the time step is halved and equations resolved. If the 
error is excessively small, the current results are accepted, stress 
and strain histories are updated but the time step for the next 
solution is increased by an arbitrary factor of 1.50. 

The following definition of error is adopted for checking equili- 
brium imbalance at the midstep or the end of the time step. For the 
i-th degree of freedom the equilibrium imbalance is given by 8S/8q ei . 

The imbalance is weighted depending upon the relative magnitudes of 
displacements. Thus, each gradient component is multiplied by the 
corresponding displacement component and the resulting work-like 
quantity is normliazed with respect to the current value of the 
potential function S. Thus the i-th component of the error vector 
is defined to be 


E. 

l 


( 9S , 

3q . q ei 


ei 


(7-1) 


The norm of the error vector is defined to be 

| ] E j | = max (E^) • (7-2) 

i 

The errors at the beginning and at the end of the time step are 
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assumed to be due solely to convergence limits of the minimization 
process. Hence the allowable error at midstep due to truncation is 
assumed to be 

INI + IWI t+At 

|| E || At — ^ (7-3) 

0 2 

At being the size of the time step. It is this measure of error at 
midstep which is required to fall within user specified limits. 
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Section 8 


APPENDIX A - FRAME ELEMENT CROSS-SECTION DETAILS 


A. 1 The Boat Cross-Section [15] 

In this section formulas for the box cross-section are given. The 
box cross-section is a thin-walled symmetric closed section composed of 
four pieces as shown in Fig. A-l. Each piece has a uniform thickness 
with the twoside pieces having the same thickness and the top and 
bottom pieces having the same thickness. 

Because of the significant difference in the torsional rigidities 
of open and closed cross-sections, use of the box cross-section to 
represent an open section by omitting pieces will lead to gross errors 
in torsional response. Thin-walled open sections can be adequately 
modeled by making a recourse to the IE section to be described later. 
The x, y, z axes of Fig. A-l are the x^, y 2 z 2 axes °f Fig* 3-3. 

The geometric parameters of the box cross section are given by 


A = 2(d 1 t 1 + d 2 t 2 ) 

y = 0 , z =0 
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yz 


0, J = n 


2 2 
2d l d 2 

V ^2 

C 1 fc 2 


(A-la, g) 


where A is the cross-sectional area, y and z are co-ordinates of the 

c c 

centroid, I , I and I are area moments of inertia, and J is the 
y yz z 

torsional constant. From symmetry the position of the shear center is 
at the centroid of the section, hence 

y s = 0, z g = 0 (A-lh,f) 
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Shear constants k and k for the cross-section are determined based on 

y z 

the maximum strain method [33], Thus 

k = ( t k - (r 

y V v y J centroid, z V v z } centroid 

y z 

Calculating the centroidal shear stresses for elastic bending the 
above formulas reduce to 
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y 2 


With these k and k , the elastic shear deflection parameters a and a 
y z y z 

of Eq. (3-34a,b) are 
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For monitoring inelastic material response, five points are spaced 
equally along each of the ,four walls of the section dividing the box 
element into 16 volume elements for the calculation of strain energy 
using Lobatto quadratures. For reporting stress results the normal 
stress, a, and the shear stress, t, are determined at the 16 points 
in Fig. A-2 and the stresses are subscripted by the point numbers 
shown in the figure to identify their location. 

For the box section, Eq. (3-44a) for the shear flow constant 
becomes 
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where the sum is over the 16 volume elements. Based on the centroidal 


strains, the measurements of the effective bending shear strains are 
Av 


Y sy L 2G (T 11 “ T 3‘ ) 


Aw g i 

Y sz = "IT = 2G (T 15 “ T 7 ) 


(A-4a,b) 


All of the equations used to treat the frame member with box section 
are now available. 


A. 2 The IE Cross-Section [15] 

In this section formulas for the IE cross-section are given. The 
IE cross-section is shown in Fig. A-3. It is made up of seven 
pieces — five flange pieces indicated by the subscripts 1, 2, 4, 6, 

7 and two web pieces indicated by the subscripts 3, 5. Each piece 
has a uniform thickness which may be different from the thickness of 
any other piece. Any combination of the five flange pieces may be 
omitted by specifying the width or thickness, or both, to be zero. 

The web sections, indicated by subscripts 3 and 5, cannot be omitted. 
The x, y, z axes of Fig. A-3 are the x^, y z ^ axes of Fig. 3-3. 

The geometric parameters of the IE cross-section are given by 
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X yz 2 A l d l + W d 3 “ ( “ A 6 d 6 + “ Ay c Z ( 


J = ^ l A - t • 

3 ii 

i=l 


(A-5a,g) 


The elastic shear center is located as follows. For loading in 

the x-z plane the normal stress is given by 

I z - I y 
z yz 

a Y M y 

I I - I y 
y z yz 

where is the bending, moment about the y-axis. From Eq. (3-24) the 
shear flow is 


q “ q n = * 


dM 


II - I 0 
y z yz 


rr f (I z - I y) tds 
2 „ z yz J ' 


where V = - y *- is the resultant shear force. This result can be used 
z dx 

to find the shear flow in each piece. 

From the definition of the shear center, we have 

d l d 2 d 6 d 7 

V z (y c- y s> = f 0 d 3 q l dS l ' f 0 d 3 q 2 ds 2 " ; 0 d 5 q 6 ds 6 + [ d 5 q 7 dS 7 


where y^ is the y co-ordinate of the shear center as shown in Fig. A-4. 
Substituting for the shear flow and integrating the following expression 


by y results 
s 


y = y - 
s c 


— i-j- {- I z [(z«3)d 3 (dJ tl -d 2 2 t 2 ) 

y z yz 
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- (y c - | VVeS, + (y c + f d 7 )d 5 d 7 t 7» 

A similar treatment for loading in the x-y plane gives an expression for 
the z co-ordinate of the shear center 


z = z - 
s c 


2 , * I y^ y c 3 d l )d 3 d l t l ^ y c + 3 d 2 )d 3 d 2 t 2 


2(1 I - I ) 
y z yz 


‘ ^c-lWeV (y c + I d 7 )d 5 d ? t 7 ] 
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- I yz [(z c + d^djt^d’t,) - (z c -d 5 )d 5 (d^ V d2t 7 )]} 

The expressions for the shear flow in each piece can he used to deter- 
mine the shear constants k and k required for the shear deflection 

y z 

calculations. Based on the maximum strain method formulas for these 


constants are 
A 


k y V ^ T y^ centroid* k 2 V ^ centroid* 
y z 


Calculation of the shear stress at the centroid based on elastic bend- 
ing in the y-z plane and substitution into the above expression yields 


for y^ positive: 
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and for y^ negative 
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Fig. A-5 Quadrature Points for IE Cross-Section. 
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The shear stress at the shear center is used to measure the effective cen- 
troidal stress, as shown in Fig. A-4 for the case of positive y . 

A similar treatment for loading in the x-z plane gives the follow- 
ing expressions for the shear constant k^: 
z c positive: 


k = 
z 


“ 2 77 {I z [A 6 +A 7 + \ *"5 ^ d 5 -z c^ ^ ^ d 5~ z c^ 


(II - I 

y z yz 5 


+ V'V/ ' y c> + A 7 ( T + ye:* + t 5 (d 5- z c )y c ]} 


z^ negative: 
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k = 
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(I 


2 ^^1^2 + 2 t 3^ d 3 +Z c^ ^ ^ d 3 +z c^ 


y z yz 3 
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+ V [A 1 ( 2“ " y c > - A 2 ( T + V " t 3 (d 3 +z c )y c 1} (A " 6c) 

With k and k determined, the elastic shear deflection parameters 
y z 

and c* z of Eqs. (3-34a,b) can be calculated. These are 
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(A-6d) 


Equations (A-5) through (A-6) represent the principal cross-section and 
member parameters used in the analysis of the IE section. They are 
based upon geometry and the elastic deflection responses. 

For monitoring inelastic material response, stress reference points 
are selected so that the end and mid-points of each of the seven cross- 
section pieces are represented. This selection divides the IE element 
into 14 volume elements for the calculation of strain energy through the 
use of Lobatto quadratures. The normal stress, O is calculated at each 
reference point as defined in Fig. A-5. The shear stress, T, is deter- 
mined only at selected reference points as shown in Fig. A-5. The shear 
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stress at the free ends of the flanges is zero. 

It is evident that the shear flow theory does not contain torsional 
components and hence cannot describe torsion of an open cross-section 
like the IE cross-section. Thus, the torsional response is assumed to 
be elastic and Eq. (6-13) is used to account for the torsional strain 
energy. 

For treating inelastic bending shear deformations, measurements of 
the effective shear strains are required. Based on the maximum strain 
method the centroidal shear stresses may be used to estimate effective 
shear strains. However, in general the centroidal stresses are a 
result of bending about both the y and z axes; hence it is necessary to 
separate the centroidal stresses into components associated with bending 
about the y-axis and with bending about the z-axis. For the box, tube 
and elliptical sections this is done easily based on geometric symmetry 
of the sections and the results are obtained in terms of centroidal 
stresses. For the non-symmetrlc IE section the analysis is more com- 
plex. Therefore a slightly different approach is taken. From Eqs . (3- 
29a, b) , the effective strains are 


V V 

Y sy ky AG ’ Y sz '^z AG 

In these expressions, the parameters, and V z are the force resultants 

of the shear stresses in the y and z directions. The forces V and V 

y z 

are determined by integrating the shear stresses over the cross-section. 
Assumption of a linear variation of shear stress between the stress 
reference points, yields the following results: 
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6 + 2 )d 3 t 3 + ( 2 + 11 + 


2 )d 5 t 5 ](A-7b) 


The shear constants and in the above equations are those deter- 
mined on the basis of elastic bending, i.e., as given by Eqs. (A-6). 

All of the equations required for the analysis of a frame member 
with IE cross-section are now available. The procedure for evaluating 
the strain energy is the same as that described for the box section. 


A. 3 The Circular Tube Cross-Section [15] 

In this section formulas for the circular tube cross-section are 
given. The circular tube cross-section is described by the mean dia- 
meter, D, and the constant wall thickness, t, as shown in Fig. A-6. The 
x, y, z axes of Fig, A-6 are the x 2 » y 2 » z 2 axes of Fig. 3-2. 

The geometric parameters of the circular tube cross-section are 
given by 



(A-8a»b) 


where A is the cross-sectional area, y and z are the co-ordinates of 

c c 

the centroid, I , I and I are area moments of inertia and J is the 
y z yz 

torsional constant. As a result of symmetry the shear center coincides 
with the centroid of the section, hence 

y g = z g = 0 (A-8c,d) 

The shear constants k and k determined on the basis of the maximum 

y z 

strain method are 
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k = (t ) „ . ,, k = “ (t ) 4 J 

y V y centro xd z V z centroid 
y z 

Calculating the centroidal shear stresses for elastic bending and 

substituting into the above formulas gives 

k =2 and k = 2 (A-9a,b) 

y z ’ 

With these k^ and k z known, the elastic shear deflection parameters 

a and a of Eqs. (3-25) are 
y z n 
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Equations (A-8) and (A-9) constitute the principal cross-section and 
member parameters used in the analysis of the circular tube cross- 
section based on elastic deflection response. 

For monitoring inelastic material response, 16 stress reference 
points spaced equally around the circumference as shown in Fig. A-6 are 
used. These points divide the circular tube section into 16 equal 
volume elements which are used for the calculation of the strain energy 
by Lobatto quadratures. These 16 reference points are used for report- 
ing the normal stress, a, and the shear stress, T, which are subscripted 
using these reference numbers as shown in Fig. A-6. 

Some of the general expressions developed for the inelastic ma- 
terial response simplify considerably for the circular tube section. 
Equation (3-40) for the shear flow may be simplified to 


= t (DG y _ i_ y 
^2L - i* L 


16 


16 


i=l 


V 


(A-10) 


where the sum is over the 16 reference points at the p-end of the 
element. Based on centroidal shear strain, measurements of the effec- 
tive shear strains caused by bending are 
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(A- 11) 


Av s ^ Aw g ^ 

Y sy = 1 T = 2G ^ T 13 _T 5' ) ’ T sz = ~L~ = 2G ^ T l" T 9^ 

All of the equations used to treat the frame member with circular tube 

cross-section are now available. 


A. 4 The Elliptical Tube Cross-Section 

In this section formulas for the elliptical tube cross-section are 
given. The elliptical tube cross-section is described by the two 
diameters 2a and 2b along the major and minor axis respectively and the 
constant wall thickness, t, as shown in Fig. A-7. The x, y, z axes of 
Fig. A-7 are the X£, y 2 » z 2 axes Fig* 3-2. 

The geometric parameters of the elliptical tube cross-section are 
A = 7T (a+b) t, y^ = - 0 


O *** 

1^ = [(3b+a)a + (3a+b) — ] 

2 

I z = ^ [(3a+b)b 2 + (3b+a) 

I =0 
yz 

T = 4(Trab) 2 t 
E(e,TT/2) 


(A-12a, f ) 


where 


e = 


AT " 2 

/a -b 


and E(e,Tr/2) is the elliptic integral 


of the second kind. 

The shear center coincides with the center as a consequence of the 
symmetry of the cross-section and hence 


y ~ z =0 
J s s 


(A-12g) 


The shear constants k and k determined on the basis of the maximum 

y z 

strain method are 


94 




Equations (A-13) through (A-14) constitute the principal cross-section 
and member parameters used in the analysis of the elliptical tube 
cross-section based on elastic deflection response. 

For monitoring inelastic material response, 16 stress reference 
points spaced equally around the circumference as shown in Fig. A-7 are 
used. These points divide the elliptical tube cross-section into 16 
equal volume elements which are used for the calculation of the strain 
energy by Lobatto quadratures. These 16 reference points are used for 
reporting the normal stress, a, and the shear stress, T, which are 
subscripted using these reference numbers as shown in Fig. A-7. 

Some of the general expressions developed for the- inelastic ma- 
terial response simplify for the elliptical tube cross-section also. 


Equation (3-40) for the shear flow may be simplified to 


t [2TTab 


16 

Ji T * ] 


4aE (e) 


(A-15) 


Where E(e) is the elliptical integral of the second kind and where the 
sum in Eq. (A-15) is over the 16 reference points at the p-end of the 
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element. Based on the centroidal shear strain, measurements of the 
effective shear strains caused by bending are 


Y 


sy 


AV s 1 AW s 1 

L = 2G ^ T 13” T 15 ) * Y sz = ~T~ = 2G ^ T l“ T 9^ 


(A- 16) 


All of the equations used to treat the frame member with elliptical tube 
cross-section are now available. 


A. 5 The Solid Rectangular Cross-Section [37] 

This section provides formulas for the solid rectangular cross- 
section which is described by its two dimensions and I >2 as shown in 
Fig. A-8. The x, y, z axes of Fig. A-8 are the x ^2* z 2 axes Fig. 
3-2. 

The geometric parameters of the rectangular cross-section are 


A = DD y = z = 0 
1 2 'c c 


D 1 D 2 D 1 D 2 

VTT- = v = 0 


(A-17a,c) 


T 1 /ri 3 X N 192 a v 1^-. tiTTb^ 
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it n=l,3... n 

Where A is the cross-sectional area, y and z the co-ordinates of the 

c c 

centroid, I , I and I are the are moments of inertia, J is the 
y z yz 

torsional rigidity constant and b is the larger and a is the shorter of 


the two dimensions and [36]. As a result of symmetry the shear 
center coincides with the centroid of the section, hence 


y = z =0 
s s 


(A-17d) 


The shear constants k and k determined on the basis of the maximum 

y z 

strain are 


k = — (r ) . , k = — (r ) . , 

y V y centroid z V z centroid 
y z 


96 




6 

8 

10 



Cj CENTROID 


w 

13 


1 * 


FIGURE A-8, 


SOLID recta: 


97 


Calculating the centroidal shear stresses for elastic bending and 
substituting into the above formulas gives 


k = k = -| (A-17e) 

y z 2 

With these k and k known, the elastic shear deflection parameters 

y z 

a and a of Eqs. (3-25) are 
y z ^ 


a 


y 


12k El 



2 * 

GAL 


a 

z 


12k El 

y_ 

2 

GAL 


(A-18a,b) 


Equations (A-17) and (A-18) constitute the principal cross-section and 
member parameters used in the analysis of the solid rectangular cross- 
section based elastic response. 

For monitoring inelastic material response, 2x2 Gauss quadrature 
points for each of the four subrectangles and 25 stress reference points 
spaced equally around the periphery as shown in Fig. A-8 are used. The 
stresses at the stress reference points are calculated by simple linear 
extrapolation of the stresses at the Gauss points. 
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